#generate function to filter out dodgy calls for GSmapper

#a single missing base
#two snps in same place

QCsnp <- function(SNPlist,sample.df=NULL){
	#SNPlist is a list of data.frames, eg hcCoveragelist2
	
	if(is.null(sample.df)) sample.df <- SNPlist[[1]]
	SNP.df <- sample.df[1,]
	for(i in 1:length(SNPlist)){
		
		tmp.df <- SNPlist[[i]]
		if(class(tmp.df)!="data.frame") next
		
		#trim out sites that appear more than once
		tmp.tab <- table(tmp.df$start.pos)
		if(any(tmp.tab>1)){
			tmp.df <- tmp.df[!tmp.df$start.pos%in%as.numeric(names(tmp.tab)[tmp.tab>1]),]
			}
		ref.nucl <- sapply(X=tmp.df$id,function(X){
			tmp <- strsplit(X,split="_")[[1]]
			tmp[length(tmp)-1]
			})
		sub.nucl <- sapply(X=tmp.df$id,function(X){
			tmp <- strsplit(X,split="_")[[1]]
			tmp[length(tmp)]
			})
		inc.site <- rep(TRUE,length(ref.nucl))
		inc.site[nchar(ref.nucl)==1&sub.nucl=="-"] <- FALSE
		inc.site[nchar(sub.nucl)==1&ref.nucl=="-"] <- FALSE
		tmp.df <- tmp.df[inc.site,]
		SNP.df <- rbind(SNP.df,tmp.df)
		}
	SNP.df <- SNP.df[-1,]
	SNP.df$cds <- unlist(SNP.df$cds)
	SNP.df$dn <- unlist(SNP.df$dn)
	SNP.df
	}
 
 
 Fst.calc <- function(p,sample.size){
	#optionally giive a vector (T/F or 1/0) giving the resistant pops
	#p is a vector of allele frequencies
	
	#do the Fst calculations (see Weir p. 147)
	if(length(sample.size)==1) sample.size <- rep(sample.size,length(sampled.mat))
	n <- sample.size
	r <- length(sample.size)
	p.hat <- sum(n*p)/sum(n)
	n.bar <- sum(n)/r
	
	sa <- sum(n*(p-p.hat)^2)/(r*n.bar) #note this is different from Weir to stop Fst>1
	
	theta <- sa/(p.hat*(1-p.hat))
	theta
	}
